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Fast Shape Estimation for Galaxies and Stars 
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ABSTRACT 

Model fitting is frequently used to determine the shape of galaxies and the point 
spread function, for examples, in weak lensing analyses or morphology studies aiming 
at probing the evolution of galaxies. However, the number of parameters in the model, 
as well as the number of objects, are often so large as to limit the use of model fitting 
for future large surveys. In this article, we propose a set of algorithms to speed up 
the fitting process. Our approach is divided into three distinctive steps: centroiding, 
ellipticity measurement, and profile fitting. We demonstrate that we can derive the 
position and ellipticity of an object analytically in the first two steps and thus leave 
only a small number of parameters to be derived through model fitting. The position, 
ellipticity, and shape parameters can then used in constructing orthonomal basis func- 
tions such as scrsiclets for better galaxy image reconstruction. We assess the efficiency 
and accuracy of the algorithms with simulated images. We have not taken into account 
the deconvolution of the point spread function, which most weak lensing analyses do. 

Key words: Galaxies: general - Methods: data analysis, statistical - Techniques: 
image processing - Gravitational lensing: shear, PSF. 



1 INTRODUCTION 

Quantifying shapes of galaxies and stars has long been one 
of the key tasks in astronomical image analyses. The shape 
of an astronomical object is one of the a few direct observ- 
ables. A lot of useful information can be inferred from the 
shapes of galaxies and stars. For example, galaxy mophology 
provides importa nt knowledge on the formation and evolu - 
tion of galaxies (|Kennicuttli 19981 : IVan de Wei et al.ll20ldl ). 
As one of the most promising tools to probe dark energy 
and dark matter in the universe, weak gravitational lensing 
relies on precision measurements of the ellipticities of the 
background galaxies and the point spread function (PSF). 

The signal in weak lensing, however, is often very weak 
and noisy, due to the large intrinsic dispersion in the shapes 
of background galaxies. For ground-based observations, the 
shapes of galaxies are further distorted by atmospheric tur- 
bulence and optical distortions, which can be described by 
the PSF, and quantified using shapes of stars. Th e first 
weak lensing algorithm was propo sed bvlKaiser et al 
(her eafter KSB) and impro ved bv lLuppino fc Kaisei 



erences therein). However, none of the algorithms satisfies 
the requirements of future surveys such as the Dark Energy 
Survey (DBS), the Canada- France-Hawaii- Telescope Legacy 
Survey (CFHTLS), and the Large Synoptic Survey Tele- 
scope (LSST). Given the enormous amount of galaxies to 
be covered by these surveys, one limiting factor of the exist- 
ing algorithms is the accuracy and efficiency with which the 
shapes and galaxies and stars are measured. 

One of the existing powerful tools that has been stud- 
ied in great detail is the decompostion of images us- 
ing basi s functions, such as th e Gaussian-Lag uerre ex- 
pans i on (IBernstein fc Jarvisll2002|'). or Shapelets (IRefregieiJ 



pans i on (Itjernstem &l JarvisllzUU2|l. or snapelets IIKetreg: 
I2OO3I : iRefregier fc Baconll2003l: iMassev fc Refregiedl2005f ) 



1995 ) 



19971 ') 



and iHoekstra et al.l (| 19981 ). Since then, a variety of algo- 
rithms were proposed and a series of data analysis chal- 
lenges have been carried out to improve the precision and re- 
duce systematic biases |Hevmans et al.l 20061: JMassev et al.l 



I2OO7I : iBridle et all |2010| . iKitching et all I2OIOI and the ref 
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serscilets ( Ngan et al.ll2009l : lAndrae et al.ll2011h . However, 
there are limitations to these methods, because that the ze- 
roth order image is often a poor match to real galaxy pro- 
files. While theoretically any intensity profile can be writ- 
ten as a weighted sum of any complete set of basis func- 
tions, in reality, the more the zeroth order resembles the 
real profile, t he fewer basis functions are needed in the 
decompostion (|Massev fc Refregieill2005l : [Boschll2010l '). As a 
first step, it is necessary to accurately and efficiently quan- 
tify the observed shape of galaxies based on noisy images, 
which can then be used as input in the construction of the 
basis functions. 

In this article, we present a set of efficient algorithms 
specifically for the task of parameter estimations for ob- 
served galaxies, including center position, ellipticity, and 
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shape parameters such as size and steepness. We then run 
numerical tests to demonstrate the effectiveness of these al- 
gorithms, and how the basis functions could help further 
improve the image reconstrcution. Our algorithms, devel- 
oped with their weak lensing applications in mind, could 
also be used for others purposes such as to measure the ap- 
proximate morphological parameters for larg e population of 
galaxies(e.g. lGadottil(2009l : IWang et "ai]|20ia ). We leave PSF 
deconvolution to later work. 

Usually, no less than six parameters are needed to model 
the light distribution of an object, including the centroid 
position {x,y), the ellipticity (31,52), the normalization Iq, 
and the profile parameter(s). High-dimensional parameter 
search is very time consuming, especially when the number 
of objects is large, and tends to be trapped at local min- 
ima. Instead of brute-force fitting, we propose to derive the 
centroid position, ellipticity, and normalization of a light 
distribution numerically and thus reduce the number of pa- 
rameters that need to be determined thr ough fitting. S imilar 
effo rt has been previousl y undertaken. iMiller et al.l (|2007l ) 
and lKitching et al.l (|2008l ) proposed a fast-fitting algorithm 
in Fourier space, which also takes into account the effects 
of PSF. Here, we will only focus on how to reproduce the 
observed images. 

Once the observed shape of galaxies (and stars for PSF 
determination) is described by smooth model(s), one may 
proceed to use the inva riants equation to arr i ve at the intrin- 
sic sh ape of galaxies (JFlusser fc SukI 119981 : iMelchior et al.l 
l2011f ). or to do t he PSF deconvolution using basis function 
transformations (JMassev fc Refregien uOOSi : IMelchior et al.l 
l2009l : lBoschll2010f ). In practice, the simple models used to de- 
scribe the light profile of stars and galaxies often lead biases 
in th e measurements l|Voigt fc Bridle! I2OICI : IMelchior et al.l 
[2OI0I). The situation can be improved by adopting more so- 
phisticated spatial models or a set of basis functions. Our 
algorithms will be usef ul in providing i nput to the construc- 
tion of basis functions (|Li fc Cuill2012r ). 

This article is organised as follows. We describe the al- 
gorithms for position and ellipticity determination in Sec- 
tions [2] and (3] and numerical tests in Section |4l We conclude 
by discussing the limitation and prospects of our algorithm 
in Section [Sj 



2 CENTROIDING 

For simplicity, we assume a Gaussian intensity profile and 
choose the origin of the coordinate system to be at the 
estimated center of the profile. With the actual center at 
{Ax, Ay), the intensity profile is given by 



I{x,y) = Ioe ^^ +e, 



(1) 



where e is the noise term. Now, with a Gaussian weight 
function, we define 



N = 



< x> 



<y> 



I{x,y)e ^"l dxdy, 


(2) 


I{x,y)e ^"ii xdxdy. 


(3) 


x2+„2 





I{x,y)e 



ydxdy. 



(4) 



Neglecting the noise terms, we have 
< a:: > 



iV 

<y> 

N 



kAx, 
kAy, 



(5) 
(6) 



where k = a'^/{a'^ -\- aj,). 

When the intensity profile has the same shape as the 
weight function, a = CTu., we have k = 0.5. If the former is 
extremely compact (i.e., o ~ 0), we have fc = 1. In general, 
the value of k is determined iteratively. For a true coefficient 
0.5 < ktrue < 1, the choice of fc = 0.75 gives |Ai+i| < ||Ai|, 
where Ai is the distance between the centers of the intensity 
profile and the weight function in the ith iteration. So, the 
calculation converges very quickly. In typical cases, conver- 
gence is achieved in less than five iterations, and the result 
can be further improved using an elliptical weight function 
with ellipticity obtained using the algorithm described in 
the next section. 



3 ELLIPTICITY MEASUREMENT 

Our approach is base d on the KSB 
shear measurement (iKaiser et al 



19951. 



algorithm for 
hereafter 



KSB95; iLuppino fc Kaiseij Il997l : iHoekstra et al.l 1 19981 ) 
Because our goal here is only to estimate the ellipticity 
of an object, there is no convolution or deconvolution 
involved. This means that we are only concerned with shear 
polarizability P^" . 

The KSB algorithm is perturbative in nature and is ap- 
plicable only when the ellipticity is small. To adapt it for 
general ellipticity measurement, we have adopted an itera- 
tive approach. At each step, we evaluate 



^g ^ 9 - gu 



(7) 



where g and gw are the ellipticities at the current and previ- 
ous steps, respectively. The difference is then added to g^ to 
obtain a better estimate of g. The iteration continues until 
the difference is less than 10^*. 

In contrast to the KSB algorithm, we use an elliptical 
weight function here, with an ellipticity also of g^. It can be 
shown (see Appendixl^ that the shear polarizability tensor 
with elliptical weight function can still be written in the 
form 



a/S 



X, 



aP 



SH 

■ eaep ■ 



(8) 



with the new X^f and ef^ given in Eqs. HA5|) and Eq. 
in Appendix |XJ respectively. 

How do we estimate Ag from ^u,? As shown in Ap- 
pendixlB] with additional terms to X and e , which we 
denote as A{X^") and A(e^-^), we have 



A5 = 5 



-{i-gl){P'"r\e° 



0), 



(9) 



where e° '^ is the observed ellipticity which is calculated by 
using our elliptical weight function, and Cu, is the ellipticity 
corresponding to g^. Eq. (|5| makes it possible to measure 
the ellipticity of an object by iterating on the ellipticity of 
the weight function unt i l it m atches that of the object. 

iBernstein fc Jarvid (|2002h proposed a similar idea but 
implemented differently. They adopted a round weight func- 
tion and try to shear the image back to be round with —gw 
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In that case, some knowledge of the multipoles of the in- 
tensity profile and t heir derivatives (which are not easy to 
compute) are needed. iHirata fc Seliald ()2003l ) also calculated 
the centroid and the ellipticity using an adptive elliptical 
weight function. Our tests show that their algorithm takes 
much more iterations to converge than ours especially for 
ellipticity because of the higher order corrections taken into 
account in our algorithm. 



4 NUMERICAL TESTS 

4.1 Galaxy Images 

To test the methods, we have created a set of galaxy images 
based on the Sersic profile. 



/(r) = Joe 



-H^V' 



(10) 



where b is chosen such that r^ is the half-light radius. There- 
fore, only three independent parameters are required to 
specify the profile. For a given galaxy, a total of seven param- 
eters are needed, including its centroid position (2 parame- 
ters) and ellipticity (2 parameters). The Sersic profile is one 
of the most popular functional forms that are used to model 
the light distribution of galaxies (e.g., iMacArthur et al.l 
I2OO3I ). 

For the simulation, we let r^ vary between 2 and 15 
pixels and n vary between 0.5 and 5. The ellipticity g is 
limited to be no greater than 0.75. For small galaxies, 2 < 
re < 10, we further require n < 0.5(re — 2) + 1 to take into 
account PSF effects, which tend to flatten the inner profile. 
We also let the centroid position vary, with respect to the 
origin, in the range of Ax, Ay G ( — 20, 20). Gaussian noise is 
added to maintain a constant signal-to-noise ratio (= 1000) 
across an image. A total of 500 512 x 512 pixels images were 
produced for the tests. 

For each galaxy, we apply the methods described in 
Sees. [2] and [3] to derive four of the seven parameters. Fur- 
thermore, the normalization factor may be derived from x^ 
minimization. 



E - 



rprcd 



(11) 



For the general case of If^^ = loff^^ , the minimization 



prcd 



dio 



leads to 



/o = 



e: 



=0 



i./rv^^^ 



ES^°'(/r")V'^? 



(12) 



(13) 



Therefore, there are only two independent parameters 
(n and re for the Sersic profile) that need to be de- 
rived from model fitting. This greatly reduces the com- 
putati on time. We used t he MINUIT minimization pro- 
gram (jjames fc Rooslll975l ) in our implementation. 

As a starting point, we estimate the half-light radius 
(rso) by locating the brightest pixel in an image and sum- 
ming up the values of neighboring pixels going outward, un- 
til the flux reaches half of the total flux. A square region 
with sides six times rso is then cut out from the 512 x 




10" 
MINUITfit:Ar/r 



Figure 1. Accuracy of ccntroiding. This shows a comparison be- 
tween our approach (fast) and the brute-force method, based the 
Sersic profile. 



512 grid. Subsequent image processing is carried out on the 
sub-image. For each galaxy, the centroid position is initially 
taken as the brightest pixel in the image. We then go through 
the iterative processes, as described in Secs.[2]and[3l to find 
the best centroid position and ellipticity. The starting val- 
ues of gi and 52 are estimated by the original KSB method. 
Iteration stops when the fractional change between two con- 
secutive steps is less than 10~^ or the number of iterations 
reaches 20. 

To assess the accuracy and efficiency of our meth- 
ods, we compare the results to those derived from the 
brute-force approach (i.e., 7-parameter fitting). Fig. [T] 
show a comparison of centroiding, Ar/re, where Ar = 
^(Aa;«=t - Ax)2 + (Ay==' - Ay)^ and {Ax,Ay) is the real 
center. Figs. [2] - H] show comparisons on Ag/g, Are/re, 
and An/n, where Ag = g''^* — g, Are = rj:^* — re, and 
An — n'^^^ — n. The superscript, "est" indicates the corre- 
sponding values derived from our algorithms and brute-force 
fitting. The accuracy of our methods is slightly worse than 
that of brute-force fitting. This can partially be attributed to 
the weight function used in our case (but not in brute-force 
fitting), as the outer region may play a non-negligible role. 
For more than 98% of the cases, Ar/re from our method 
ranges between 0.01% and 1%. 

Our estimation of the ellipticity deviates from the true 
ellipticity by less than 1% in most cases. Like any moments- 
based ellipticity estimation algorithms, our estimation of p is 
biased when g is large, due to the limited size of the image 
stamps used (side length Grso), resulting in moments cal- 
culations not converging completely. However, this bias is 
significantly reduced if we use our measured shape parame- 
ters to construct sersiclets basis functions. The reproduced 
image is extroplated to a bigger size (side length lOrso) so 
that subsequent moments calculations are not effected by 
the cut-off boundary. Then we run our algorithms again on 
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Figure 2. Accuracy of cUipticity measurements. This siiows 
a comparison between our approach (fast) and the brute-force 
method, based the Sersic profile. The dots are obtained using cl- 
lipticity algorithms decribed in Sec.|3l and crosses are results after 
using the nine lowest order sersiclets basis functions constructed 
on a larger pixel grid. The side length of the larger pixel grid here 
is chosen to be ten times the half-light radius. 
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e e 

Figure 3. Accuracy of half-light radius determination. This 
shows a comparison between our approach (fast) and the brute- 
force method, based the Sersic profile. 

the sersiclets modellj The crosses in Fig. [5] shows the Ag/g 
estimated using the sersiclets model with the nine lowest or- 
der basis functions build on a larger pixel grid. If we continue 
this process, each time we construct a new sersiclets model, 
the data points in Fig. [2] continue to move toward the diago- 
nal line until it reaches an unbiased state. Fig. [S] shows how 

^ Such operations should bo paid more attention to when the 
cUipticities of the outter and inner regions differ significantly. 




-0.1 

MINUITfit :An/n 



0.1 



Figure 4. Accuracy of Sersic index estimation. This shows a com- 
parison between our approach (fast) and the brute-force method, 
based the Sersic profile. 



/^g/g varies with the true eUipticity for brute-force fitting, 
our algorithms, and our algorithms combined with sersiclets 
modelling. The effectiveness of our approach in removing the 
bias at large eUipticity is clear. We also note that the wider 
spread of points in Fig. [S] (c) than in (a) is because that 
the model used in brute-force fitting is the same as what 
we used to simulate the images, in which case the brute- 
force fitting is superior over any other methods in terms of 
accuracy. But for real galaxies that are typically more com- 
plex than a simple Sersic model (e.g., with bulge, bar, spiral 
arms, etc.), the basis functions will be more adva ntageous 
fe.g.. lMassev fc Refregieill2005l : lAndrae et al.ll201ll ). 

Because the centroid position and eUipticity are both 
fixed in the subsequent two-parameter fitting, the accuracy 
on re and n is worse than that of the brute-force approach, 
but is still well within 10% in most cases, as shown in Fig. [3] 
and Fig. |4l Most of the the outliers have small re and big n. 
Their images are dominated just by several central pixels. 
This leads to relatively large uncertainty on the centrioding 
and then biases re and n bacause the strong degeneracy 
between these two parameters. 

Going beyond individual parameters, we assess the 
quality of the fit to the overall light profile of a galaxy. Fig. [6] 
shows a comparison of the reduced x^ values. When Afpixoi 
is large, the ratio x^/-^pixei is expected to be around unity 
if the estimation is unbiased. We found that the fit becomes 
progressively worse, as the Sersic profile steepens (i.e., large 
n) and the size decreases. When the half-light radius ap- 
proaches the size of a pixel, small error in centroiding may 
contribute significantly to the x^ value. As one would have 
expected from the results in Fig. [2] constructing sersiclets 
model using the estimated parameters reduces x^, as shown 
by the crosses in Fig. [6] 

As for efficiency, on a single-core 2.2 GHz Intel CPU, it 
takes about 0.1 seconds to carry out centroiding and eUiptic- 
ity measurement on one galaxy image and about 0.5 seconds 
to derive the remaining two parameters in the Sersic model 
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Figure 5. The dependence of ellipticity estimation accuracy on 
tlic size of ellipticity for (a) brute-force fitting, (b) our algorithms, 
and (c) our algorithms combined with sersiclets modelling. 
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Figure 6. A comparison of the reduced x^ values between our 
method (fast) and the brute-force method, based on the Sersic 
profile. The dots are obtained using algorithms described in 
Sees. [2] and |3] and crosses are results after using the nine low- 
est order sersiclets basis functions. 



from profile fitting. In comparison, it takes about 6 seconds 
to process one galaxy image in the brute-force approach. 



4.2 Star Images 

We carried out similar tests with simulated star images, 
which are relevant to PSF modelling. Here, instead of us- 
ing fitting programs like MINUIT, we developed faster nu- 
merical methods also to derive shape parameters from the 
moments of observed light profiles. Two kind of PSF profiles 
are considered here. 



4-2.1 Gaussian profile 
For the Gaussian profile, 

I{r) = loe'^^ , 
the zeroth and second order moments are 



ro 



r2 = 



I{r)W{r)27vrdr, 



I{r)W{r)2Tvr-^dr, 



where the weight function is 





W{r) = e 2" 


If we take a^ - 


= a, we have 




ro — TTCT lo 




r2 — TTCr Jo 


Therefore, 










ro 



/o 



(14) 

(15) 
(16) 

(17) 

(18) 
(19) 

(20) 
(21) 



For the simulation, the half-light radius (rso) is allowed 
to vary from 1.5 to 2.5 pixels. The ellipticity g is limited to no 
greater than 0.25. We also let the centroid position vary, with 
respect to the origin, in the range of Ax, Ay £ (—20,20). 
Gaussian noise is added to maintain a constant signal-to- 
noise ratio (= 300) across an image. A total of 500 512 x 
512 images are produced for the tests. 

As for galaxies, we carry out centroiding and ellipticity 
determination with the methods described in Sees. [2] and 
O We then use Eqs. PUI) and (PT]) to derive a and Jo from 
the zeroth- and second-order moments. Setting a^ = cr, we 
recalculate the moments. The iteration continues until the 
fractional change in a between two consecutive steps is less 
than 10^ or the number of iterations reaches 20. 

Fig. [7] shows a comparison between our approach and 
brute-force fitting in determining the half-light radius. In 
practice, the integrations of Eqs. (fT5|l and (fT6|l were done 
numerically. And this will bias the result because of pixe- 
lation, which has a prominent effect on small images. How- 
ever, deviation is still within a few percent. The bias in the 
overall size can be roughly estimated as Sr^o/ < rfg > ^ 
0.001. This value is acceptabl e for weak lensing statistics 
(jPaulin-Henriksson et al.ll2009l ). The overall quality of the 
shape estimation is shown in Fig. [8] 



6 G. L. Li, B. Xm& W. Cm (2012) 

0.04r 

0.03 

0.02- 



0.01 - 



-0.01 



-°l°^ 




02 -0.01 



0.01 0.02 0.03 0.04 



MINUITfit:Ar^„/r^„ 

OU DU 



Figure 7. Half-light radius. This shows a comparison between our 
approach (fast) and the brute-force method, based the Gaussian 
profile. 
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Figure 9. A comparison of the determination of half-light radius 
between our method (fast) and the brute-force method. This is 
based on the Moffat profile. 
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Figure 8. A comparison of the reduced x^ values between our 
method (fast) and the brute-force method, based on the Gaussian 
profile. 



4.2.2 Moffat profile 
For the MofTat profile 



I[r) = Jo 1 + 



Td 



in addition to ro and r2, we define 

ri= I I{r)W(r)2'n:r^dr, 



(22) 



(23) 



where tlie weight function W{r) is also of the MofTat shape, 

-3 



W{r) = 1 + 



Tdn 



If we take rd-w = rd, we have 



r2 = /o 
Vi = 7o 



4 



(/3 + l)(^ + 2)' 
/3(^ + l)(^ + 2)- 



Therefore, 









Io = 



r-o(/J + 2) 



(24) 

(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



As for the Gaussian profile, we have made 500 512 x 
512 images. In the simulation, /3 varies between 3 and 5. 
Similarly, we compute the centroid position, ellipticity, and 
rd from the moments by iterating on rdw. The results are 
shown in Figs, (gj - (fTT)) . 

Again, the biases in rso and /? are attributable to pix- 
elation effects when we do the numerically integration of 
moments. The half-light radius rso is generally accurate to 
about 3-4% and the power-law index j3 to less than 20%. 
As before, the bias in the overall size is roughly Sr^o/ < 
'^50 > ^ 0.001. Overall, the shape estimation is worse than 
in the case of Gaussian profiles (comparing Fig. [11] with 
Fig. |8]). We found that the points with x^ > 1-5 are dom- 
inated by realizations with r^ < 2. This is because that 
the pixelation effect bias the moments two much. As same 
reason for Sersic model, the Moffat profile also suffers from 
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Figure 10. A comparison of the determination of Moffat power- 
law index between our method (fast) and the brutc-forcc method. 
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Figure 11. A comparison of the reduced x^ values between our 
method (fast) and the brute-force method, based on the Moffat 
profile 



rd-P degeneracy. Therefore the problem here is more tough 
than that in Gaussian profile. In practice, for the data taken 
under similar weather conditions, we suggest to use several 
bright stars to derive r^ and P by using fitting process and 
then fix /3 at the average value for other (fainter) stars. This 
could be a way out to break the degeneracy and to decrease 
the pixelation effect. 

Based on the shape estimation algorithms described and 
tested above for the Gaussian and Moffat profiles, we have 
constructed two set of basis functions which are named gaus- 
sianlets and moffatlets, whose z eroth order profi les are Gaus- 
sian and Moffat, respectively (jLi fc Cuill2012l ). The corre- 



sponding algorithms have been tested in the GRavitational 
lEnsing Accuracy Tes ting 2010 (GREATIO) Star Challenge 
(jKitchine: et al.l 120101 ) . Our gaussianlets worked very well 
with accuracy of a{e) ~ a( R^)/R^ ~ 0.1%, wher e e is the 
ellipticity and R is the size (jKitching et al.ll2012l ). But the 
moffatlets didn't do very well especially when the PSF size 
is small. 



5 DISCUSSION 

In this work, we have adapted the KSB algorithm for ellip- 
ticity determination. The new algorithm allows the use of 
elliptical weight functions and is applicable to highly ellipti- 
cal shapes. This, combined with centroiding, makes it pos- 
sible to eliminate four of the parameters from brute-force 
fitting. Consequently, less time is needed in searching for 
the remaining (shape) parameters. Overall, the efficiency is 
improved roughly by an order of magnitude. 

We have tested our algorithms with simulated images 
that are representative of stars and galaxies. In general, 
our centroiding algorithm is accurate to well within 1%, al- 
though it is slightly worse than the brute-force fitting. For 
galaxy images, which are made with the Sersic profile, our 
approach is relatively worse than brute-force fitting in de- 
termining the half-light radius and the Sersic index. This is 
mainly due to the fact that the two parameters are tightly 
coupled, and then the fitting results will be biased by any 
small offset in centroiding. Overall, we can find very good 
fit to simulated galaxy profiles with our algorithms (as mea- 
sured by the reduced x^)- The accuracy of parameter esti- 
mations is seen to be further improved when the algorithms 
are run on sersiclets models with nine lowest order basis 
functions which are constructed with the estimated param- 
eters as input. Given that galaxies are not perfectly Sersic- 
like, our algorithms, combined with sersiclets modelling, has 
unique advantages over other shape estimation techniques. 

For star images, we have adopted profiles that are of- 
ten used for PSF modelling, Gaussian and Moffat. In this 
case, we show that we can also derive the shape parameters 
directly from the moments of the light distribution, which 
improve the efficiency even further. However, for small im- 
ages, pixelation leads to significant biases. Fortunately, the 
bias in the overall size is Sr^o/ < rfo > < 0.001, which is still 
acceptable for weak lensing statistics. We participated in the 
GREATIO star challenge. The results show that our meth- 
ods performed quite well with Gaussian profiles but not as 
satisfactorily with Moffat profiles, particularly for data sets 
with stars of small radii. This tells us that the effects of pix- 
elation and rd-P degeneracy must be taken into account for 
Moffat profile. 

Like other algorithms, we also assume axisymmetric 
light distributions with constant ellipticity. In practice, this 
is known not to be the case for either galaxies or PSFs. 
In fact, the morphology of a galaxy or PSF can be quite 
complex. Therefore the use of simple spatial models leads 
to biases. The accuracy is expected to improve with more 
complicated spatial models. One possible option is to adopt 
appropriate basis functions, such as sersiclets, moffatlets and 
gaussianlets, to reconstruct the light distribution. The ap- 
proach described here shows its strong capability in provid- 
ing input parameters for the basis functions. In spite of the 
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deficiencies, our algorithms are very fast, which is impor- 
tant for future surveys, as the data volume is expected to 
increase drastically. 
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APPENDIX A: THE KSB ALGORITHM WITH AN ELLIPTICAL WEIGHT FUNCTION 

Lets W{R^) be an elliptical weight function, which is assumed to be converted by shear polarization guj, 

R" = [(1 - g-.,if + gl2]0l + [(1 + g^if + gl2]0l - ^9^26162. (Al) 

Eq. (B5) in KSB95 now takes on the form 



Zlmii — 



dOi 



= W{s,i6j6^ + s^ie^Ora + 5mie,e,) + w'e,e,era^^. (A2) 

OOl 

Compared to the cases with circular weight function, the only change is that 26i is replaced with !„ ' in the terms that 



2 



contain W' , where the prime denotes differentiation with respect to R , i.e., 

01 ^ i^^ = Ae^ - 2ff»2e2, (A3) 

02^^^^^Be2~2g^2Oi, (A4) 

Z (7(72 

where A — 1 + g^ — 2g^i and _B = 1 + g™ + 2gu,i. X and e are now given by 

^SH _i [ a^fif,.( 2W6'' + 2W\el~el)(Ael-Bel) 2W' {el -el)[{A + 3)9^62 ~2g^2e^] \ 

- T j ^^ >\ AW'e^e2{A6l - Bel) "iwe^ + AW'6^e2[{A + B)e^e2 - 2g,2e^] ) ' ^'^^> 

and 

'" ^ < :: ) + 1 /"■'"/<«)»■'''" ( ,A+ BilsZ-t-o' ) ■ <'^"'* 

APPENDIX B: THE ITERATIVE PROCESS FOR ELLIPTICITY DETERMINATION 

The KSB method is applicable only when the ellipticity of an object is small, i.e., (^,ij in Eq. (B7) in KSB95 is small, so 
that the Taylor expansion is appropriate. However, galaxies are often highly elliptical. This appendix shows a derivation of 
Agr = g — g^^ where g is the intrinsic ellipticity of an object and g^ is the estimated ellipticity. 
Let /3 be the coordinates in the virtual plane with ellipticity (;,„. (3 and 6 are related by 

f3 = CO. (Bl) 

where 

^ ^ f 1 - "ti - (gu^ili + g^2'y2) -12 + gw2ii - gn,i'y2 \ ,g2) 

V -72 - gw2ii + 3uji72 1 + 71 - (511.171 + gn,2'y2) J ' 



-^I^' --T^' 



and 



The surface brightness is conserved during the transformation defined in Eq. (IB2[) . 

/„(/3) = /(0). (B4) 

Let ip — C~^ — 1, it follows that 

/„(/3) = /(C7-^/3) = /(/3 + W). (B5) 

Following the original KSB formalism, 

Qij ~ Qij ~ WlmZlmij, (B6) 

where Zimij is defined in Eq. (B4) in KSB95. In their Eqs. (3.2) and (3.3), the difference in ellipticity e between the real image 
with reduced shear g and the virtual image with g^ is easily related to SQij — Qij — QYj- We then find that, in addition to 
the terms defined in Eqs. (|A5|) and (|A6|) in Appendix[Xl we now have the following extra terms on X^^ and e'^^: 

^ ' T J ^^ '\ -51P3 - 52P4 giP4 -~g2P3 J ' ^ ' 

and 

A(e-)^i/.^./(.)(;f^-f^«), (B8) 
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where 



Pi = AWiel - el) + 2W'{el - el){Ael + Bel - '^gie^ei), (B9) 

Pa = 'iwe^ei + 2W'{el - el)[{A - B)eie2 + 2g2{ef - el)], (bio) 

P3 = 8weie2 + 4W'eie2{Aej + Bel - 4520102), (bii) 

P4 = -2W{el - el) + 4W'eie2[{A - P)6»i6»2 + 252(0? - el)], (B12) 

Ps = 4VK0^ + 2VK'0^(A6'? + Bel - 4520102), (B13) 

Pe = 2W'e^[{A - P)0i02 + 252(0? - el)]. (B14) 
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